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1. Motivation and objectives 

The problem of transition from laminar to turbulent flow in free shear layers is of 
great practical interest. Many natural and industrial situations involve a free shear 
flow, and it is very crucial to understand the mechanisms governing the process of 
transition to turbulence in order to predict and, if possible, control the evolution 
of the flow. Much past work related to the instability of the mixing layer has dealt 
exclusively with Newtonian fluids. The interaction between the large scale two- 
dimensional structures, namely the roll-up and the pairing, and vortex stretching 
has been explored extensively, leading to a better understanding of the mechanisms 
of transition to turbulence. 

In a wide range of applications, the reduction of hydrodynamic drag is of great im- 
portance for reasons of economy and performance. Delaying or accelerating laminar- 
to- turbulent transition of a free shear layer has many obvious advantages. To delay 
transition to turbulence as far downstream as possible allows a gain in energy effi- 
ciency in some industrial processes. Accelerating the transition, on the other hand, 
can be of interest in processes where high mixing is desired. 

Among the modern techniques for drag reduction are riblets, large-scale eddy 
breakup devices, and flow suction. Another avenue for drag reduction is the use of 
polymer additives in the flow. The study of drag reduction obtained by the addition 
of small amounts of high polymers has been an active area of research for the last 
three decades (see Sellin and Moses 1989). 

This work is a continuing effort to advance our understanding of the effects of 
polymer additives on the structures of the mixing layer. In anticipation of full non- 
linear simulations of the non-Newtonian mixing layer, we examined in a first stage 
the linear stability of the non-Newtonian mixing layer (Azaiez and Homsy 1992). 
The results of this study show that, for a fluid described by the Oldroyd-B model, 
viscoelasticity reduces the instability of the inviscid mixing layer in a special limit 
where the ratio (^) is of order 1 where We is the Weissenberg number, a measure 
of the elasticity of the flow, and Re is the Reynolds number. 

In the present study, we will pursue this project with numerical simulations of 
the non-Newtonian mixing layer. Our primary objective is to determine the effects 
of viscoelasticity on the roll-up structure. We will also examine the origin of the 
numerical instabilities usually encountered in the simulations of non-Newtonian 
fluids. 
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2. Accomplishments 

2.1 Problem Definition 

We used a vorticity-streamfunction formulation for Cauchy’s momentum equa- 
tion. This equation is closed through evolution equations relating the stress tensor 
to the shear rate tensor. In all the subsequent analysis, the stress tensor is written 
as the sum of two stresses: 

t~t s + t p ( 1 ) 

In Eq. (1), the first term corresponds to the contribution of the Newtonian solvent 
and is proportional to the shear rate tensor: 

r * = *7*7 (2) 

where rj 9 is the solvent viscosity. The second term represents the polymeric contri- 
bution: 

t p = T] p a ( 3 ) 

where rj p is the polymeric contribution to the shear viscosity and a is the tensor 
related to the polymeric stress. Equations for a(x,y,t) define a rheological model 

and must be given to close the equations. We let k = — — — = — be the ratio 

V* + Vp V 

of viscosities. We fix k = 0.5 in all the simulations. 

The first step in studying the non-Newtonian behavior of a given flow is the 
choice of the appropriate constitutive equation required to describe the physics of 
the non-Newtonian fluid. In the present study, we will focus on two rheological 
models known to describe dilute polymer solutions: The popular Oldroyd-B model 
and the FENE model. 

We are initially interested in the roll-up of the shear layer into spanwise vortices. 
The initial base state for the numerical simulations is given by a hyperbolic tangent 
velocity profile on which we superposed a cos perturbation on the stream function 
with the maximum at the centerline y — 0. The most unstable wave number, 
a — 0.44, has been used in all the simulations. The extension of the domain in the 

27T 

streamwise direction is set to — 6 with 6 the momentum thickness of the mixing 

a 

layer. The vertical extension has been heuristically fixed equal to 86. The Reynolds 

number is defined as Re = — - where u 0 = (U\ — £M/2. U\ and Uo are the 

rj 

free-stream velocities. 

The vorticity equation and the rheological equations lead to a system of partial 
differential equations that we solve by a pseudo-spectral method based on the Hart- 
ley transform. We use a fourth order Runge-Kutta method to advance the equations 
for the Hartley coefficients in time. 

2.2 The Oldroyd-B model 

The Oldroyd-B model describes well the behavior of polymeric liquids composed 
of a low concentration of high molecular weight polymer in a very viscous Newtonian 
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solvent at moderate shear rates. These fluids have come to be known as Boger 
fluids. Data on well-characterized Boger fluids are available (Boger 1977; Mackay 
& Boger 1987), allowing a comparison of theoretical predictions with experimental 
results. The Oldroyd-B rheological equations can be derived from a molecular 
model in which the polymer molecule is idealized as a Hookean spring connecting 
two Brownian beads (Bird et al. 1987). 

As above, the stress tensor t can be written as the sum of two stresses: 

T = r* + T p = rjlicj + (1 - «)a] 

The tensor a satisfies the upper convected Maxwell equation: 

a,j + = 7 «> ( 4 ) 


where: 

£^LL = ^ +i ;.Va-V{r- L a-a-Vt7 (5) 

St dt 

is the upper- convected derivative of a, and A is the polymer relaxation time. Using 
u 0 and 8 as the reference velocity and the reference length, respectively, Eq. (4) is 

\u Q 

characterized by one dimensionless group: the Weissenberg number, We — -j-, a 
dimensionless measure of the elasticity of the flow. The Oldroyd-B model contains 
both the upper- convected Maxwell fluid (rj, = 0 k = 0) and the Newtonian fluid 
(rj p = 0 k = 1 ). 

For small values of the Weissenberg number, We, the flow does not show any 
changes from the Newtonian case. Obviously, viscoelastic terms are not large enough 
to affect the inertial mechanisms. 

As we increase the Weissenberg number, the code starts to develop instabilities 
and ceases to converge before the completion of the roll-up. Similar problems re- 
lated to the complexity of the coupling between viscoelastic and inertia effects have 
been reported by Joseph et al. (1985) and Crochet and Delvaux (1990). Numerical 
instabilities in the simulations of viscoelastic fluids have in fact been the subject 
of many studies. These studies have in general focused on the Maxwell model, a 
special case of the Oldroyd-B model, and showed that the absence of the diffusive 
term in the vorticity equation may be the reason for the blow-up of the numer- 
ical schemes. Other studies attribute the breakup of the numerical codes to the 
reaching of a supercritical regime where the vorticity equation becomes hyperbolic 
(Joseph (1990)). This regime is reached whenever the viscoelastic Mach number 
M = (We Rt) 1 ! 2 exceeds one. (Crochet and Delvaux (1990)) 

In the present study, we shall show that the divergence of codes using the Oldroyd- 
B model is attributed to a relaxation term in the stress equations which plays a 
major role in the development and growth of the instability. Writing Eq. (5) in 
component form, we have: 


da\2 

dt 
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da 11 

dt 


. 1 d 2 xj> . . d d , n ,d 2 xy tp d 2 t/> 

” ^~Wt +2 lhd^ au ~^ u !h + v d^ an + 2 ^~wT +ai2 W ^ 
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l* 1 relaxation coeff. 
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Damping Term 


Convective Term 


da 
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dt 
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■We dxdy 
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d_ 
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, dhp 

We +au dx 2 > 
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2 nd Relaxation Coeff. 

v 

Damping Term 


Convective Term 


We analyzed the evolution of the different terms in the stress equations for various 
values of the Weissenberg number. We find that the relaxation coefficient behaves 
differently as we change the Weissenberg number. 
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FIGURE 1. Time evolution of the first relaxation coefficient. 


The “relaxation coefficients” are spatially dependent and coupled to the vorticity 
transport equation. Fig. 1 shows the evolution of the minimum and the maximum 
of the first relaxation coefficient for three different values of We. The blow-up of 
the code is always associated with positive values of this coefficient. (It is worth 
mentioning that for another rheological model, the Jeffrey’s corotational model, we 
were able to conduct successful simulations for very large value of We. In fact the 

relaxation coefficients for this model are all negative and equal to — — 

W e 
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The vorticity. 
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FIGURE 2. Vorticity and 1st Relaxation Coefficient Contours for the Oldroyd-B 
Model. Re=100,We=5. 

The contour plots presented in Fig. 2 show that the regions where the first relax- 
ation coefficient becomes positive correspond to the regions of subsequent numerical 
instability. In these regions, the an component of the stress tensor starts to build 
up, leading to a blow-up of the code. The same results are obtained for the second 

relaxation coefficient. As a conclusion, it resorts from this analysis that in regions 

dxi Ov 

of the flow where the parameter M = We|( — ~ ^)l 6 ets larger than one > the 

stresses start growing causing the code to blow up. The parameter M involves the 
elasticity of the polymer and the local normal stress in the flow. 

2.8 The F.E.N.E model 

The Oldroyd-B model gives a steady state elongational viscosity that goes to 
infinity at a finite elongational rate. This unlikely behavior results because the 
Hookean dumbbell model permits infinite extension. In order to avoid this unreal- 
istic behavior, a Warner law is used instead of the Hook law leading to the FENE 
model. This model is described by Mackay and Petrie (1989) and the rheological 
equation is: 

a,,Z + We 6 -^-=y, 1 + ~(l + a, i We) (6) 

where Z = 1 + — (1 + a n ), n = 2 for a two dimensional flow. The FENE model 

is characterized with a third parameter b , related to the nature of the spring used 
to model the macromolecule. When b — > oo this rheological model reduces to the 

Oldroyd-B model. In the numerical solution of Eq. (6), —— was derived from 

the transport equation for an . The finite extension properties of the FENE model 
allows full nonlinear simulation of roll-up. 
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FIGURE 3. Evolution of the energy of the perturbation. 

The evolution of the energy in the fundamental is presented in Fig. 3. The 
region of linear growth shows that the non-Newtonian mixing layer described by 
the FENE model is more unstable than the Newtonian mixing layer. As we shall 
see later, the roll-up occurs faster in the non-Newtonian case and the structure of 
the mixing layer changes from the Newtonian case. This higher instability of the 
non-Newtonian mixing layer agrees with our linear stability analysis (not discussed 
here). 

Fig. 5 shows contour plots of the vorticity for Re = 100, We = 50, and b = 5. 
Notice the changes in the structures of the roll-up as compared to the Newtonian 
case (Fig. 4). The average value of the vorticity does not differ much from the 
Newtonian case, but its distribution has clearly changed. Intense accumulation of 
vorticity tends to occur in the braids and continue into a spiral form within the 
core of the vortex. The intense gradients of vorticity persist for a longer time in the 
core, leading to a faster and more intense roll-up. 

3. Conclusion and future plans 

The simulation of the non-Newtonian mixing layer was very instructive and seems 
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FIGURE 4. Vorticity contours for the Newtonian fluid. Re— 100. 
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FIGURE 5. Vorticity contours for the FENE model. Re— 100, We— 50, and b— 5. 

to be very promising for reaching a better understanding of the mechanisms of insta- 
bility of the mixing layer and for controlling them through polymer additives. For 
the Oldroyd-B rheological model, the numerical code ceases to converge whenever 
one of the relaxation coefficients in the stress equations becomes positive. These co- 
efficients involve the elasticity of the polymer and the extensional stress in the flow. 
We propose to pursue this analysis in order to understand the physical meaning of 
these coefficients in terms of the coupling between viscoelastic and inertia effects. 
The numerical simulations using the FENE model lead to very interesting results. 
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For the parameters studied, free shear flows described by this rheological model are 
more unstable than their Newtonian counterparts and we obtained a fairly good 
description of the effects of viscoelasticity on the roll-up. The process of roll-up for 
these flows occurs quickly and induces intense gradients of vorticity in the braids 
and the core of the flow. The next step of the present study is to study the effects 
of viscoelasticity on pairing. We will then proceed to examine the effects on vortex 
stretching and the streamwise vortices of the three dimensional mixing layer. 
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